

	
*****************************************************************		
* be careful not to accidentally overwrite this 	
preserve 
		clear 
		set obs 1 
		gen var1=. 
		
		save "$datapath/intermediate/collect_robustness_sigma.dta", replace		
restore 
*****************************************************************		




	
*****************************************************************		
** Define the program that'll maximize librarian utility

capture program drop find 

program define find
		tempvar Dgp 
		tempvar Dge 
		tempvar sDg 
		tempvar Qp
		tempvar Qe
		tempvar b
		tempvar c
		tempvar d
		tempvar e
		tempvar f
		tempvar ff
		tempvar g
		tempvar con
		tempvar ne
		tempvar np
		tempvar n2

     if "`3'"!=""{
        loc sse = "`3'"
        qui cap drop `sse'
        }
     else {
        tempvar sse
        }
		

			gen double  `Dgp' = (abs(`1'[1,1])*exp($dpvar/(1-$sbvar)))^((1-$sbvar)/(1-$stvar)) 
			gen double `n2' = ($bvar - abs(`1'[1,1])*$ppvar)/$pevar 
			gen double `Dge'  = (abs(`n2')*exp($devar/(1-$sbvar)))^((1-$sbvar)/(1-$stvar))
			gen double `sDg'  = `Dgp' + `Dge' 
			gen double `Qp'  = $mvar*`Dgp'*(`sDg'^(-$stvar))/(1+`sDg'^(1-$stvar))
			gen double `Qe'  = $mvar*`Dge'*(`sDg'^(-$stvar))/(1+`sDg'^(1-$stvar))
			gen double `np' = abs(`1'[1,1])
			gen double `ne' = abs(`n2')
			
			gen double `f' = ($tvarp*`Qp' + $tvare*`Qe')  - 1000000*abs( $bvar - $pevar*abs(`n2') - $ppvar*abs(`1'[1,1]) )
			gen double `con' = $bvar - $pevar*abs(`n2') - $ppvar*abs(`1'[1,1])
			su `f' `con'
			su `np' `ne'
			
			egen double  `sse' = sum(`f')
		
			sca `2' = `sse'
		
		
	end
*****************************************************************		

	


use $datapath/clean/main_data.dta, clear 


    
*****************************************************************		
** Create the sigma and delta variables from demand estimation

	reshape long Q N, i(popu_lsa year lno ) j(type) s
	
	gen double M = 50*popu_lsa
	gen double s = Q/(N*M) 
	
	drop if s==.
	keep if s>0 & s<1 
	
	gen double de = type=="e"
	egen double QQ = sum(Q), by(lno year)
	
	gen double ls = ln(s) - ln((M - QQ)/M)
 
	gen double lstop = ln(Q/QQ)
	egen double NN = sum(N), by(lno year)

	gen double lstop_inst = ln(N/NN)
	gen double lsbottom = ln(1/N)
	gen double lsoverall = ln(1/NN)
	
	gen bsbottom = .
	gen bstop = .
	
tempfile loopdata
save `loopdata'
	

foreach sigmatop in 0.1 0.4 0.8 {
	foreach sigmabottom in 0.5 0.937 {
		
		if `sigmatop' < `sigmabottom' {
 	
			use `loopdata', clear
			
			replace bsbottom =`sigmabottom'
			replace bstop = `sigmatop'
			
			
			gen double del = ls - bstop*lstop - bsbottom*lsbottom   
		*****************************************************************		

			keep if year==2018 
	
	
		*****************************************************************		
		** prepare variables (derivatives) needed to calculate the thetas
	
			gen double Dg = (N*exp(del/(1-bsbottom)))^((1-bsbottom)/(1-bstop))
			egen double sDg=sum(Dg), by(lno year)
			gen double shat = (1/N)*Dg*(sDg^(-bstop))/(1+sDg^(1-bstop))
			
			gen double N1p = N + .1*(1-de)
			gen double N1e = N + .1*de 
			
			gen double Dg_1p = (N1p*exp(del/(1-bsbottom)))^((1-bsbottom)/(1-bstop))
			egen double sDg_1p=sum(Dg_1p), by(lno year)
			gen double shat_1p = (1/N1p)*Dg_1p*(sDg_1p^(-bstop))/(1+sDg_1p^(1-bstop))

			gen double Dg_1e = (N1e*exp(del/(1-bsbottom)))^((1-bsbottom)/(1-bstop))
			egen double sDg_1e=sum(Dg_1e), by(lno year)
			gen double shat_1e = (1/N1e)*Dg_1e*(sDg_1e^(-bstop))/(1+sDg_1e^(1-bstop))

	
		* note: dqp is the derivs of the two q's wrt p 
			gen double dqp = M*(N1p*shat_1p - N*shat)/.1 
			gen double dqe = M*(N1e*shat_1e - N*shat)/.1 

		
		** Set prices 
			** materials costs calculated as mean expenditure/mean holdings from table 1 
				
			gen pricep = 0.91
			gen pricee = 0.61
			
			gen pp = pricep
			gen pe = pricee
		

*****************************************************************		

	

	
*****************************************************************		
** Calculate thetas based on sigmas (and relative prices)

	preserve 
	
		keep lno type dqp dqe pp pe 
		reshape wide dqp dqe, i(lno ) j(type) string
		rename  dqpp dqp_dnp  
		rename  dqee dqe_dne  
		rename  dqep dqp_dne  
		rename  dqpe dqe_dnp  
		
		gen double det = dqp_dnp*dqe_dne - dqp_dne*dqe_dnp 
		
		gen double thetap = (dqe_dne/det)*pp - (dqe_dnp/det)*pe 
		gen double thetae = -(dqp_dne/det)*pp + (dqp_dnp/det)*pe 
		
		keep lno theta*
		reshape long theta, i(lno) j(type) string 
		drop if theta==. 
		egen double nn=count(theta), by(lno)
		keep if nn==2 
		save theta.dta , replace 
		
	restore 

	merge 1:1 lno type using theta.dta  
	keep if _merge==3
*****************************************************************		
	
	
*****************************************************************		
** Turn it into a wide dataset 	

	gen price = pp
	replace price= pe if type=="e"

	keep N theta bstop bsbottom price M del lno type Q 
	
	reshape wide N theta price del Q, i(lno bstop bsbottom M) j(type) string 
	
	gen ratio = thetae/thetap 
	su ratio, de 
	
	gen double budget = Ne*pricee + Np*pricep
	
*****************************************************************		
	
	

	

*****************************************************************		
*****************************************************************		
** Find optimal mNp and mNe for each library
*****************************************************************		
*****************************************************************		


* randomize the order 
	gsort lno 
	set seed 86
	drawnorm rando 
	gsort -rando 

local mmm =_N  

forvalues j= 1 (1) `mmm' {	


************************************		
** baseline and double Pe


	foreach f in 10 20 { 
	
	preserve 

	    replace pricee =pricee*(`f'/10)
	
		keep if _n==`j'
	
		gl nlist = "beta:A"	
		gl dpvar = "delp"
		gl devar = "dele"  
		gl bvar = "budget"  
		gl stvar = "bstop"  
		gl sbvar = "bsbottom"
		gl ppvar = "pricep"  
		gl pevar = "pricee"  
		gl tvarp = "thetap"  
		gl tvare = "thetae"  
		gl mvar = "M"  
		gl np0var = "Np"  

		su Np 
		local p=r(mean)
		
		su Ne
		local e=r(mean)

				 
		 mat b0 = ( `p' )
		 mat colnames b0 = $nlist
		 
		 amoeba find b0 yout bols . . 1E-10 
		 
*****************************************************************		

		 
		 
*****************************************************************		
** Calculate holdings, quantities and CS based on optimal holdings

		 
		gen mNp = bols[1,1]
		gen mNe = (budget - pricep*mNp)/pricee 
		
		su Np Ne mNp mNe 
		gen factor = `f'
			gen double  Dgp = (mNp*exp(delp/(1-bsbottom)))^((1-bsbottom)/(1-bstop)) 
			gen double Dge  = (mNe*exp(dele/(1-bsbottom)))^((1-bsbottom)/(1-bstop))
			gen double sDg = Dgp + Dge 
			gen double mQp  = M*Dgp*(sDg^(-bstop))/(1+sDg^(1-bstop))
			gen double mQe  = M*Dge*(sDg^(-bstop))/(1+sDg^(1-bstop))
			gen double CS  = M*log(1+sDg^(1-bstop))
		
*****************************************************************		
	

*****************************************************************		
** save the variables we care about

		keep Np Ne mNp mNe lno factor Qe Qp  mQe mQp CS lno M dele delp thetae thetap bstop bsbottom 

		gen raise ="Pe"
		
		sleep 200
		capture append using "$datapath/intermediate/collect_robustness_sigma.dta" 
		capture save  "$datapath/intermediate/collect_robustness_sigma.dta", replace  

	restore 
		
	}

**************************************	
	
	

************************************		
** double Pp
	
	preserve 

		replace pricep = 2*pricep
	
		keep if _n==`j'
	
		gl nlist = "beta:A"	
		gl dpvar = "delp"
		gl devar = "dele"  
		gl bvar = "budget"  
		gl stvar = "bstop"  
		gl sbvar = "bsbottom"
		gl ppvar = "pricep"  
		gl pevar = "pricee"  
		gl tvarp = "thetap"  
		gl tvare = "thetae"  
		gl mvar = "M"  
		gl np0var = "Np"  

		su Np 
		local p=r(mean)
		
		su Ne
		local e=r(mean)

		 mat b0 = ( `p' )
		 mat colnames b0 = $nlist
		 
		 amoeba find b0 yout bols . . 1E-10 
		 
*****************************************************************		

		 
		 
*****************************************************************		
** Calculate quantities and CS based on optimal holdings

		 
		gen mNp = bols[1,1]
		gen mNe = (budget - pricep*mNp)/pricee 
		
		su Np Ne mNp mNe 
		gen factor = 20
			gen double  Dgp = (mNp*exp(delp/(1-bsbottom)))^((1-bsbottom)/(1-bstop)) 
			gen double Dge  = (mNe*exp(dele/(1-bsbottom)))^((1-bsbottom)/(1-bstop))
			gen double sDg = Dgp + Dge 
			gen double mQp  = M*Dgp*(sDg^(-bstop))/(1+sDg^(1-bstop))
			gen double mQe  = M*Dge*(sDg^(-bstop))/(1+sDg^(1-bstop))
			gen double CS  = M*log(1+sDg^(1-bstop))
		
*****************************************************************		
	

*****************************************************************		
** save the variables we care about

		keep Np Ne mNp mNe lno factor Qe Qp  mQe mQp CS lno  M dele delp thetae thetap bstop bsbottom
		
		gen raise ="Pp"
		
		
		
		sleep 200		
		capture append using "$datapath/intermediate/collect_robustness_sigma.dta" 
		capture save   "$datapath/intermediate/collect_robustness_sigma.dta", replace  

	restore 
		
			}	

	

		}
	}
}
